HS project

Trajectory from HF-SCs using TInGa

This script is used to infer trajectory from HF-SCs, using TInGa.

library(dplyr)
library(patchwork)
library(ggplot2)

.libPaths()

## [1] "/usr/local/lib/R/library"

Preparation

In this section, we set the global settings of the analysis. We will store data there :

save_name = "non_matrix"
out_dir = "."

We load the sample information :

sample_info = readRDS(paste0(out_dir, "/../../1_metadata/hs_hd_sample_info.rds"))
project_names_oi = sample_info$project_name

graphics::pie(rep(1, nrow(sample_info)),
              col = sample_info$color,
              labels = sample_info$project_name)

This are the settings for trajectory inference :

traj_dimred = "harmony"      # to infer trajectory on
traj_2d = "harmony_dm"       # just of visualization
seed = 1337L
traj_max_dims = 40           # to infer trajectory on
name2D = "harmony_18_tsne"   # just of visualization

We load the Seurat object :

sobj = readRDS(paste0(out_dir, "/", save_name, "_sobj.rds"))
sobj

## An object of class Seurat 
## 17837 features across 5599 samples within 1 assay 
## Active assay: RNA (17837 features, 2000 variable features)
##  8 dimensional reductions calculated: RNA_pca, RNA_pca_18_tsne, RNA_pca_18_umap, harmony, harmony_18_umap, harmony_18_tsne, harmony_dm, harmony_dm_5_umap

Set root

We set the root manually :

root_cell_id = "2022_14_GAACACTAGTGCCTCG-1"
sobj$is_root = colnames(sobj) == root_cell_id
sobj$cell_name = colnames(sobj)
root_plot = aquarius::plot_label_dimplot(sobj, reduction = traj_2d,
                                         col_by = "is_root", col_color = c("gray80", "red"),
                                         label_by = "cell_name", label_val = root_cell_id) +
  ggplot2::ggtitle("Root cell") +
  ggplot2::theme(plot.title = element_text(hjust = 0.5, face = "bold"),
                 aspect.ratio = 1) +
  Seurat::NoLegend() + Seurat::NoAxes()

sobj$is_root = NULL
sobj$cell_name = NULL

root_plot

Trajectory inference

Now, we perform the trajectory inference :

set.seed(seed)
my_traj = aquarius::traj_tinga(sobj,
                               seed = seed,
                               expression_assay = "RNA",
                               expression_slot = "data",
                               count_assay = "RNA",
                               count_slot = "counts",
                               dimred_name = traj_dimred,
                               dimred_max_dim = traj_max_dims,
                               root_cell_id = root_cell_id,
                               tinga_parameters = list(max_nodes = 3))

## Add pseudotime
sobj$pseudotime = my_traj$pseudotime

(Time to run : 9.31 s)

Visualization

On the diffusion map

We visualize pseudotime :

Seurat::FeaturePlot(sobj, features = "pseudotime",
                    reduction = traj_2d,
                    cols = viridis::viridis(n = 100)) +
  ggplot2::lims(x = range(sobj@reductions[[traj_2d]]@cell.embeddings[, 1]),
                y = range(sobj@reductions[[traj_2d]]@cell.embeddings[, 2])) +
  ggplot2::labs(title = "Pseudotime") +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_text(hjust = 0.5)) +
  Seurat::NoAxes()

We visualize the trajectory with arrows :

dynplot::plot_dimred(trajectory = my_traj,
                     label_milestones = FALSE,
                     size_milestones = 0,
                     plot_trajectory = TRUE,
                     dimred = sobj[[traj_2d]]@cell.embeddings,
                     color_cells = 'pseudotime', color_trajectory = "none")

We color cells according to cluster ID :

dynplot::plot_dimred(trajectory = my_traj,
                     label_milestones = FALSE,
                     plot_trajectory = TRUE,
                     dimred = sobj[[traj_2d]]@cell.embeddings,
                     grouping = sobj$seurat_clusters,
                     size_milestones = 0,
                     color_trajectory = "none")

On the tSNE

We visualize pseudotime :

Seurat::FeaturePlot(sobj, features = "pseudotime",
                    reduction = name2D,
                    cols = viridis::viridis(n = 100)) +
  ggplot2::lims(x = range(sobj@reductions[[name2D]]@cell.embeddings[, 1]),
                y = range(sobj@reductions[[name2D]]@cell.embeddings[, 2])) +
  ggplot2::labs(title = "Pseudotime") +
  ggplot2::theme(aspect.ratio = 1,
                 plot.title = element_text(hjust = 0.5)) +
  Seurat::NoAxes()

We visualize the trajectory with arrows :

dynplot::plot_dimred(trajectory = my_traj,
                     label_milestones = FALSE,
                     size_milestones = 0,
                     plot_trajectory = TRUE,
                     dimred = sobj[[name2D]]@cell.embeddings,
                     color_cells = 'pseudotime', color_trajectory = "none")

We color cells according to cluster ID :

dynplot::plot_dimred(trajectory = my_traj,
                     label_milestones = FALSE,
                     plot_trajectory = TRUE,
                     dimred = sobj[[name2D]]@cell.embeddings,
                     grouping = sobj$seurat_clusters,
                     size_milestones = 0,
                     color_trajectory = "none")

Save

We save the Seurat object :

saveRDS(sobj, file = paste0(out_dir, "/", save_name, "_sobj_traj_tinga.rds"))

We save the trajectory object :

saveRDS(my_traj, file = paste0(out_dir, "/", save_name, "_my_traj_tinga.rds"))

R Session

show

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
## 
## locale:
## [1] C
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] dynutils_1.0.5  dynwrap_1.2.1   purrr_0.3.4     tidyr_1.1.4    
## [5] ggplot2_3.3.5   patchwork_1.1.2 dplyr_1.0.7    
## 
## loaded via a namespace (and not attached):
##   [1] softImpute_1.4              graphlayouts_0.7.0         
##   [3] pbapply_1.4-2               lattice_0.20-41            
##   [5] haven_2.3.1                 dyndimred_1.0.3            
##   [7] vctrs_0.3.8                 usethis_2.0.1              
##   [9] blob_1.2.1                  survival_3.2-13            
##  [11] prodlim_2019.11.13          later_1.3.0                
##  [13] DBI_1.1.1                   R.utils_2.11.0             
##  [15] SingleCellExperiment_1.8.0  rappdirs_0.3.3             
##  [17] uwot_0.1.8                  dqrng_0.2.1                
##  [19] gng_0.1.0                   jpeg_0.1-8.1               
##  [21] zlibbioc_1.32.0             pspline_1.0-18             
##  [23] pcaMethods_1.78.0           mvtnorm_1.1-1              
##  [25] htmlwidgets_1.5.4           GlobalOptions_0.1.2        
##  [27] future_1.22.1               UpSetR_1.4.0               
##  [29] laeken_0.5.2                leiden_0.3.3               
##  [31] clustree_0.4.3              lmds_0.1.0                 
##  [33] parallel_3.6.3              scater_1.14.6              
##  [35] irlba_2.3.3                 DEoptimR_1.0-9             
##  [37] tidygraph_1.1.2             Rcpp_1.0.9                 
##  [39] readr_2.0.2                 KernSmooth_2.23-17         
##  [41] carrier_0.1.0               promises_1.1.0             
##  [43] gdata_2.18.0                DelayedArray_0.12.3        
##  [45] limma_3.42.2                pkgload_1.2.2              
##  [47] graph_1.64.0                RcppParallel_5.1.4         
##  [49] Hmisc_4.4-0                 fs_1.5.2                   
##  [51] RSpectra_0.16-0             fastmatch_1.1-0            
##  [53] ranger_0.12.1               digest_0.6.25              
##  [55] png_0.1-7                   sctransform_0.2.1          
##  [57] cowplot_1.0.0               DOSE_3.12.0                
##  [59] here_1.0.1                  TInGa_0.0.0.9000           
##  [61] dynplot_1.1.0               ggraph_2.0.3               
##  [63] pkgconfig_2.0.3             GO.db_3.10.0               
##  [65] DelayedMatrixStats_1.8.0    gower_0.2.1                
##  [67] ggbeeswarm_0.6.0            iterators_1.0.12           
##  [69] DropletUtils_1.6.1          reticulate_1.26            
##  [71] clusterProfiler_3.14.3      SummarizedExperiment_1.16.1
##  [73] circlize_0.4.15             beeswarm_0.4.0             
##  [75] GetoptLong_1.0.5            xfun_0.35                  
##  [77] bslib_0.3.1                 zoo_1.8-10                 
##  [79] tidyselect_1.1.0            GA_3.2                     
##  [81] reshape2_1.4.4              ica_1.0-2                  
##  [83] pcaPP_1.9-73                viridisLite_0.3.0          
##  [85] rtracklayer_1.46.0          rlang_1.0.2                
##  [87] hexbin_1.28.1               jquerylib_0.1.4            
##  [89] dyneval_0.9.9               glue_1.4.2                 
##  [91] waldo_0.3.1                 RColorBrewer_1.1-2         
##  [93] matrixStats_0.56.0          stringr_1.4.0              
##  [95] lava_1.6.7                  europepmc_0.3              
##  [97] DESeq2_1.26.0               recipes_0.1.17             
##  [99] labeling_0.3                httpuv_1.5.2               
## [101] class_7.3-17                BiocNeighbors_1.4.2        
## [103] DO.db_2.9                   annotate_1.64.0            
## [105] jsonlite_1.7.2              XVector_0.26.0             
## [107] bit_4.0.4                   mime_0.9                   
## [109] aquarius_0.1.5              Rsamtools_2.2.3            
## [111] gridExtra_2.3               gplots_3.0.3               
## [113] stringi_1.4.6               processx_3.5.2             
## [115] gsl_2.1-6                   bitops_1.0-6               
## [117] cli_3.0.1                   batchelor_1.2.4            
## [119] RSQLite_2.2.0               randomForest_4.6-14        
## [121] data.table_1.14.2           rstudioapi_0.13            
## [123] org.Mm.eg.db_3.10.0         GenomicAlignments_1.22.1   
## [125] nlme_3.1-147                qvalue_2.18.0              
## [127] scran_1.14.6                locfit_1.5-9.4             
## [129] scDblFinder_1.1.8           listenv_0.8.0              
## [131] ggthemes_4.2.4              gridGraphics_0.5-0         
## [133] R.oo_1.24.0                 dbplyr_1.4.4               
## [135] BiocGenerics_0.32.0         TTR_0.24.2                 
## [137] readxl_1.3.1                lifecycle_1.0.1            
## [139] timeDate_3043.102           ggpattern_0.3.1            
## [141] munsell_0.5.0               cellranger_1.1.0           
## [143] R.methodsS3_1.8.1           proxyC_0.1.5               
## [145] visNetwork_2.0.9            caTools_1.18.0             
## [147] codetools_0.2-16            Biobase_2.46.0             
## [149] GenomeInfoDb_1.22.1         vipor_0.4.5                
## [151] lmtest_0.9-38               msigdbr_7.5.1              
## [153] htmlTable_1.13.3            triebeard_0.3.0            
## [155] lsei_1.2-0                  xtable_1.8-4               
## [157] ROCR_1.0-7                  BiocManager_1.30.10        
## [159] scatterplot3d_0.3-41        abind_1.4-5                
## [161] farver_2.0.3                parallelly_1.28.1          
## [163] RANN_2.6.1                  askpass_1.1                
## [165] GenomicRanges_1.38.0        RcppAnnoy_0.0.16           
## [167] tibble_3.1.5                ggdendro_0.1-20            
## [169] cluster_2.1.0               future.apply_1.5.0         
## [171] Seurat_3.1.5                dendextend_1.15.1          
## [173] Matrix_1.3-2                ellipsis_0.3.2             
## [175] prettyunits_1.1.1           lubridate_1.7.9            
## [177] ggridges_0.5.2              igraph_1.2.5               
## [179] RcppEigen_0.3.3.7.0         fgsea_1.12.0               
## [181] remotes_2.4.2               scBFA_1.0.0                
## [183] destiny_3.0.1               VIM_6.1.1                  
## [185] testthat_3.1.0              htmltools_0.5.2            
## [187] BiocFileCache_1.10.2        yaml_2.2.1                 
## [189] utf8_1.1.4                  plotly_4.9.2.1             
## [191] XML_3.99-0.3                ModelMetrics_1.2.2.2       
## [193] e1071_1.7-3                 foreign_0.8-76             
## [195] withr_2.5.0                 fitdistrplus_1.0-14        
## [197] BiocParallel_1.20.1         xgboost_1.4.1.1            
## [199] bit64_4.0.5                 foreach_1.5.0              
## [201] robustbase_0.93-9           Biostrings_2.54.0          
## [203] GOSemSim_2.13.1             data.tree_1.0.0            
## [205] rsvd_1.0.3                  memoise_2.0.0              
## [207] evaluate_0.18               forcats_0.5.0              
## [209] rio_0.5.16                  geneplotter_1.64.0         
## [211] tzdb_0.1.2                  caret_6.0-86               
## [213] ps_1.6.0                    DiagrammeR_1.0.6.1         
## [215] curl_4.3                    fdrtool_1.2.15             
## [217] fansi_0.4.1                 highr_0.8                  
## [219] urltools_1.7.3              xts_0.12.1                 
## [221] GSEABase_1.48.0             acepack_1.4.1              
## [223] edgeR_3.28.1                checkmate_2.0.0            
## [225] scds_1.2.0                  cachem_1.0.6               
## [227] randomForestSRC_2.12.1      desc_1.4.1                 
## [229] npsurv_0.4-0                babelgene_22.3             
## [231] rjson_0.2.20                openxlsx_4.1.5             
## [233] ggrepel_0.9.1               clue_0.3-60                
## [235] rprojroot_2.0.2             stabledist_0.7-1           
## [237] tools_3.6.3                 sass_0.4.0                 
## [239] nichenetr_1.1.1             magrittr_2.0.1             
## [241] RCurl_1.98-1.2              proxy_0.4-24               
## [243] car_3.0-11                  ape_5.3                    
## [245] ggplotify_0.0.5             xml2_1.3.2                 
## [247] httr_1.4.2                  assertthat_0.2.1           
## [249] rmarkdown_2.18              boot_1.3-25                
## [251] globals_0.14.0              R6_2.4.1                   
## [253] Rhdf5lib_1.8.0              nnet_7.3-14                
## [255] RcppHNSW_0.2.0              progress_1.2.2             
## [257] genefilter_1.68.0           statmod_1.4.34             
## [259] gtools_3.8.2                shape_1.4.6                
## [261] HDF5Array_1.14.4            BiocSingular_1.2.2         
## [263] rhdf5_2.30.1                splines_3.6.3              
## [265] AUCell_1.8.0                carData_3.0-4              
## [267] colorspace_1.4-1            generics_0.1.0             
## [269] stats4_3.6.3                base64enc_0.1-3            
## [271] dynfeature_1.0.0            smoother_1.1               
## [273] gridtext_0.1.1              pillar_1.6.3               
## [275] tweenr_1.0.1                sp_1.4-1                   
## [277] ggplot.multistats_1.0.0     rvcheck_0.1.8              
## [279] GenomeInfoDbData_1.2.2      plyr_1.8.6                 
## [281] gtable_0.3.0                zip_2.2.0                  
## [283] knitr_1.41                  ComplexHeatmap_2.14.0      
## [285] latticeExtra_0.6-29         biomaRt_2.42.1             
## [287] IRanges_2.20.2              fastmap_1.1.0              
## [289] ADGofTest_0.3               copula_1.0-0               
## [291] doParallel_1.0.15           AnnotationDbi_1.48.0       
## [293] vcd_1.4-8                   babelwhale_1.0.1           
## [295] openssl_1.4.1               scales_1.1.1               
## [297] backports_1.2.1             S4Vectors_0.24.4           
## [299] ipred_0.9-12                enrichplot_1.6.1           
## [301] hms_1.1.1                   ggforce_0.3.1              
## [303] Rtsne_0.15                  shiny_1.7.1                
## [305] numDeriv_2016.8-1.1         polyclip_1.10-0            
## [307] grid_3.6.3                  lazyeval_0.2.2             
## [309] Formula_1.2-3               tsne_0.1-3                 
## [311] crayon_1.3.4                MASS_7.3-54                
## [313] pROC_1.16.2                 viridis_0.5.1              
## [315] dynparam_1.0.0              rpart_4.1-15               
## [317] zinbwave_1.8.0              compiler_3.6.3             
## [319] ggtext_0.1.0